Nonequilibrium effects in tunnel Josephson junctions 
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We study nonequilibrium effects in current transport through voltage biased tunnel junction with long diffu- 
sive superconducting leads at low applied voltage, eV <C 2A, and finite temperatures. Due to a small value of 
the Josephson frequency, the quasiparticle spectrum adiabatically follows the time evolution of the supercon- 
ducting phase difference, which results in the formation of oscillating bound states in the vicinity of the tunnel 
junction (Andreev band). The quasiparticles trapped by the Andreev band generate higher even harmonics of 
the Josephson ac current, and also, in the presence of inelastic scattering, a nonequilibrium dc current, which 
may considerably exceed the dc quasiparticle current given by the tunnel model. The distribution of travelling 
quasiparticles also deviates from the equilibrium due to the spectrum oscillations, which results in an additional 
contribution to the dc current, proportional to \/V. 

PACS numbers: 74.45,+c, 74.40.+k, 74.25.Fy, 74.50.+r. 



I. INTRODUCTION 



Tunnel current transport in superconducting junctions is a 
classical topic of interest and research. 1 Theoretical descrip- 
tion of the superconducting tunneling given by the tunnel 
model 2 - 3 is based on the assumption about equilibrium in the 
junction electrodes. The resulting equations for the tunnel cur- 
rent are expressed through the nonperturbed density of states 
(DOS) and equilibrium distribution function. Such an ap- 
proach includes only single-particle tunnelling processes, and 
it is sufficient for describing current-voltage characteristics at 
large applied voltage eV > 2A, and within the subgap volt- 
age region eV < 2A at nonzero temperature. At very small 
temperature, multiparticle tunneling processes must be taken 
into account^ However, these processes are exponentially 
weak at small voltage, and if the temperature is not particu- 
larly small, eV <C T < 2A, only tunneling of thermally excited 
quasiparticles plays significant role. It is clear, however, that 
quasiparticle tunnelling generates nonequilibrium distribution 
in the electrodes, and that the effect is enhanced in diffusive 
junctions because of the scattering by impurities. Tradition- 
ally, this effect is considered to be negligibly small since it is 
of a higher order in the junction transparency. 

How small is the nonequilibrium effect actually? In volt- 
age biased superconducting junctions, the phase of the order 
parameter has different time dependencies in different elec- 
trodes, and the interference of the order parameters induced 
by the tunnelling leads to the time oscillations of the DOS. 
The character of these oscillations can be qualitatively under- 
stood from comparison with the well studied case of ballis- 
tic tunnel junctions. 6,7,8 In the latter, the Andreev levels are 
formed in the vicinity of the tunnel barrier; their energies lie 
within the energy gap in the bulk electrodes and oscillate in 
time following evolution of the superconducting phase differ- 
ence. Similarly, one may expect that in diffusive junctions the 
time oscillation of the DOS will have the form of a "breath- 
ing" potential well localized in the vicinity of the tunnel bar- 
rier ("Andreev band"), which periodically, with the Josephson 
period, traps and releases quasiparticles. Thus the quasipar- 
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FIG. 1: Possible experimental realizations of tunnel junctions (a,b) 
and the theoretical model (c): the tunnel barrier attached to bulk su- 
perconducting electrodes by long (L >• £q) diffusive superconducting 
leads. 

tides spend larger time within the junction area compared to 
the travelling time, which should generate strongly nonequi- 
librium quasiparticle distribution and, as a result, an appre- 
ciable change of the tunnelling current. The effect should be 
most pronounced for a planar junction geometry sketched in 
Fig-Da>b), when the tunnel junction is connected to bulk elec- 
trodes via diffusive superconducting wires whose length ex- 
ceeds the size of the Andreev band. Otherwise, as in diffusive 
point contacts, the proximity of large equilibrium reservoirs 
will suppress the DOS oscillations and hence formation of the 
Andreev band because of rapid spreading out of the current. 

In this paper we demonstrate that the nonequilibrium quasi- 
particle distribution generated by nonstationary process of for- 
mation of the Andreev band during the tunneling considerably 
modifies the dc tunnel current at small applied voltage. We 
show that the nonequilibrium effect results in a time depen- 
dent amplitudes of all the three current components given by 
the tunnel model* 2 * 3 - 

I(t) =h{t) sin <j>+I 2 (t)+l3 (f)cos0, §=2eVt (1) 

(Josephson current, quasiparticle current, and the interference 
current, respectively). All the amplitudes have nonharmonic 
time dependence, and generally all of them contribute to the 
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dc current. However, at small voltage, the major contribu- 
tion, nonanalytic injunction transparency and voltage, comes 
from the sine term (Josephson current). More precisely, this 
additional dc current results from a non-adiabatic component 
of the distribution function associated either with inelastic re- 
laxation of quasiparticles, or with the quasiparticle diffusion 
away from the tunnel barrier. The additional dc current may 
considerably exceed the magnitude given by the tunnel model, 
and moreover, leads to a non-monotonous net dc current. 

Such effect has close qualitative similarity to the nonequi- 
librium effects in transparent weak links [superconductor-nor- 
mal metal-superconductor (SNS) junctions and superconduct- 
ing bridges] studied earlier. 9,10 There are however the differ- 
ences: In superconducting bridges, the potential well appears 
due to suppression of the order parameter by the current con- 
centration (departing effect), 9 while in the tunnel junctions 
this effect does not play any essential role. In long SNS junc- 
tions, the potential well is pre-prepared by the proximity ef- 
fect, and it exists in the absence of the transport current; as a 
result, the oscillations develop at small energy, 10 while in the 
tunnel junctions the oscillations occur at large energies close 
to the gap edge in the electrodes. 

Organization of the paper is the following. In Sec. II we 
formulate the theoretical model and microscopic equations 
describing current transport through the tunnel junction with 
diffusive leads. In Sec. Ill we introduce an adiabatic approach 
for solving the time-dependent problem in the limit of low 
applied voltage. This approach is applied to calculation of the 
spectral characteristics of the junction (Sec. IV) and the quasi- 
particle distribution function (Sec. V). In Sec. VI we calculate 
the ac and dc components of the net Josephson current, and 
then discuss and summarize the results in Sec. VII and VIII. 



H. MODEL AND BASIC EQUATIONS 

The model of tunnel junction we are going to study is de- 
picted in Fig. 0c) and consists of a tunnel barrier with the 
transparency T attached to bulk superconducting electrodes 
via two superconducting leads (— L < x < and < x < L). 
We will consider diffusive limit, in which the elastic scat- 
tering length £ is much smaller than the coherence length 
<jjo = \J D /2 A, where D is the diffusion coefficient (we as- 
sume h = ks = 1). The length L of the leads is assumed to 
be much larger than <^o (long junction), and their width will 
be supposed to be much smaller than the Josephson penetra- 
tion depth which implies homogeneity of the current along the 
junction. Similar model was considered in Refs. fill andfl2l in 
study of the dc Josephson effect in tunnel structures. 

Under these conditions, the microscopic calculation of the 
electric current 7(f) requires solution of the one-dimensional 
diffusive equations of nonequilibrium superconductivity 13 
(see also a review 14 ) for the 4x4 matrix two-time Keldysh- 
Green's function G(x,t\,t 2 ) in the leads, 

[6,6] = iVdJ, J = God x G, G 2 = 8(ti-t 2 ), (2) 
ft = [ia z d ti - e<p+A(t{)]8(h-t 2 ), 



where A = e" 7 ^ 7<7 V A, <p is the electric potential, A and <j) are 
the modulus and the phase of the order parameter, respec- 
tively, (7; are the Pauli matrices, and 

<5=({ J*), G K =fo?-}of. (3 ) 

In Eq. 0, g RA are the 2x2 Nambu matrix retarded and ad- 
vanced Green's functions, and / = f + a z f- is the matrix dis- 
tribution function (we use 'check' for 4x4 and 'hat' for 2x2 
matrices). The multiplication procedure in Eqs. (0-0 is de- 
fined as the time convolution, 

(AoB)(t u t 2 ) = J A(h,t)B(t,t 2 )dt. 

In Eqs. 0, we neglect the inelastic collision term, which will 
be taken into account later, in Sec. VI. 

The boundary conditions for the function G and the matrix 
current J at the left (x = — 0) and the right (x = +0) sides of 
the tunnel junction are given by relation^ 

J-o = J+o = (2g N Ry 1 [G_oAo] , (4) 

where R is the junction resistance and gpj is the normal con- 
ductivity of the leads per unit length. The electric current is 
related to the Keldysh component of the matrix current / as 

I(t) = (7Cg N /4e)Tra z f K (x,t,t), 

and thus it can be expressed through the boundary value Jq, 

I(t) = ^Tro z [G^a,G +0 \ K (t,t). (5) 

Equations can be decomposed into the diffusion equa- 
tions for the Green's functions, 

[H,g]= iVdJ, f=go d x g, g 2 = 8(h-t 2 ), (6) 

and the equation for the Keldysh component G K , 

[ft, G K ) = Wd x f K , J K = f o d x G K + G K o d x g A . (7) 

The boundary conditions for the functions g and G K at the 
tunnel barrier follow from Eq. 

fo = (W/So)[g- ,g +0 ], (8a) 
Jj = (W/^o)[G- 0l G +0 ] K . (8b) 
In Eqs. 0, the transparency parameter W is defined as 

w = r{£d)/2R - (3<W^)r » r, 

where R(%o) — ^,o/gN is the normal resistance of the lead per 
length ^ has been shown in Ref s . ITU andlO that this quan- 
tity rather than the barrier transparency T plays the role of 
a true transparency parameter in the theory of diffusive tun- 
nel junctions (see also discussion in Sec. VII). In this paper, 
we will consider the limit W <C 1, which corresponds to the 
conventional tunneling concept. In this case, according to 
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Eqs. (|8), the gradients of all functions are small. Within the 
tunnel model, which assumes W to be the smallest parameter 
in the theory, these gradients are neglected, and the functions g 
and / are assumed to be local-equilibrium within the leads. In 
our consideration, we will lift this assumption and suppose the 
local-equilibrium form of these functions only within the bulk 
electrodes (reservoirs). Attributing the reference point for the 
phase, (j) = 0, to the left electrode, x — —L, these functions in 
the right electrode, x = L, are given by relations, 

g(E,t) = a z u s {E + a z eV) + ie ia ^a y v s (E), (9a) 
(u s ,v s )^ = JE ^ w - 2 ,E RA =E±iO, 

f(E) = n eq (E + a z eV), n eq {E) = tanh(£/2r), (9b) 

written in terms of the mixed Wigner representation A (E,t) of 
the two-time functions, 

A(h,t 2 ) = J + ~^- e - E ^A(E,t). (10) 

In Eq. (II Ot . the variable E has the meaning of the quasiparticle 
energy, and t = {t\ +?2)/2 is a real time. Similar equations, 
with = and V = 0, apply to the left electrode, x = — L. 

Because of the small value of the tunneling parameter W 
one can neglect the variation of the electric potential and the 
superconducting phase along the leads, and assume the volt- 
age V and the phase difference <p = 2eVt between the reser- 
voirs to be directly applied to the tunnel barrier. Following 
this argument, one can also neglect the variation of the charge 
imbalance function /_ proportional to a small electric field 
(~ eVW) penetrating the superconducting leads. Furthermore, 
the small value of the superfluid momentum in the supercon- 
ducting leads|ii*i£ p s ~ W/£o, enables us to neglect a small 
effect of the suppression of the energy gapjl& [~ (p s ^o) 4 ^ 3 ~ 
W 4 / 3 ], by the superfluid momentum. In such approximation, 
the coefficients in the equation (|2ji within the left lead, x < 0, 
are time-independent functions, similar to the value of G at 
the left electrode. At x > 0, using the gauge transformation, 

G(x > 0,h,t 2 ) = S f (t i )G(x > 0, 1 1 , t 2 )S{t 2 ), of the function G 
with a unitary operator S(t) = exp[/cr,0(f)/2]jli we exclude 
the time-dependent phase and the electric potential from the 

equations for the function G and the boundary conditions 
at x = L, which then become similar to the equations for G(x) 
at x < and the boundary conditions at x = —L, This results 

in the symmetry relation G(x) = G(—x), which allows us to 
replace the function G + q in the boundary condition (0J and in 
the expression (0 for the electric current by the inverse gauge 
transformation of the function Go = G_o, 

G+o Go = S(fi)(5+oS + (h) = S{h)G sHt 2 ). (11) 

As the result, the problem is reduced to the solution of a static 
equation for the function G(x) within the left lead with the 
time-dependent boundary condition (|4j at the tunnel barrier. 
Similar approach is used in the theory of ballistic point con- 
tacts, where the Josephson coupling is described by an effec- 
tive time-dependent matching condition for the gauge-trans- 
formed Bogolyubov-de Gennes equations in the leads. 



In a general nonstatic case, the function G(t\,t 2 ) consists 
of a set of harmonics G(E„,E m ), E n = E +neV, which are 
coupled to each other through a complicated set of recursive 
equations following from Eqs. (0 (see discussion in Ref.fl8i). 
The problem essentially simplifies if the distance eV between 
the harmonics is much smaller than the smallest scale 8E of 
variations in the quasiparticle spectrum, 

eV < SE. (12) 

The magnitude of 8E will be indicated below, see Eq. 

III. LOW VOLTAGE LIMIT 

In the limit of low applied voltages, the evolution of the 
quasiparticle spectrum and distribution function can be de- 
scribed within the adiabatic approximation using expansion 
over small Josephson frequency. In the static case, eV — > 0, 
the function G depends only on the time difference t\ — 1 2 , 
and the Wigner transformation JlQi reduces Eqs. (|6j to the 
standard Usadel equations^ 

Ev — Au= (iD/2)d x (ud x v — vd x u) , (13a) 

u 2 -v 2 =l, (13b) 

for the scalar components of the Green's function 

g(x,E) = G-u + iOyV. (14) 

The functions u and v determine the spectral characteristics of 
the system; in particular, the quantity N(x,E) = (u R — ur)/2 
is the DOS normalized over its value Nf in the normal state. 
In what follows, we will express the advanced Green's func- 
tions through the retarded ones, (u, v) A = — («, v) R *, using the 
general relation g^ = —a z g RJ( G z , and omit the superscript R, 
assuming all Green's functions to be retarded. 

At small applied voltages, we proceed to the Wigner repre- 
sentation d 1 Oft of the two-time functions and expand the time 
convolutions in Eqs. (|6j-(|8j to first order in eV, 

(AoB){E,t)^AB+(i/2)(A,B), 

where (A,B) = 3eA d t B — d t A 3eB denotes the Poisson brack- 
ets in the energy-time space. Within such approximation, the 
Green's function g(x,E,t) holds the matrix structure (114-1 . and 
the gauge-transformed functions g Q and f [see Eq. il 1H . 
which enter the boundary conditions l|8}, read 

1 = o z uo(E + a z eV,t) + ie 2kJ -~ eVt o y v {E,t), (15a) 

f =fo(E + a z eV,t). (15b) 

The expression for the electric current obtained from 
Eqs. (|5} and (1151 consists of the three terms, 



7(f) 


= Ii sin 0+72+73 cos <p , 


(16a) 


7i(f)sin0 


= ^J°°JEI s (E,t)ME,t), 


(16b) 


h(t) 


= £ / dEN 2 (E,t)d E f Q (E,t), 
R Jo 


(16c) 


h{t) 


= £ rdEM 2 (E,t)d E f (E,t), 
R Jo 


(16d) 
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FIG. 2: Energy dependencies of the DOS for the transparency pa- 
rameter W = 0.05: (a) at the tunnel barrier for different values of the 
superconducting phase; (b) at = % for different distances from the 
barrier. 

where fy(E,t) is the boundary value of the distribution func- 
tion f(x,E,t), and I s = -Imv'osin^, Nq = N(x = 0,E) = 
Re wo, and Mq = Re vo are the spectral densities of the different 
current components. Such structure of the current is similar to 
the result of the tunnel model: 3 indeed, when the current spec- 
tral densities and distribution function approach nonperturbed 
equilibrium form, the first term in Eq. dl6a> describes the ac 
Josephson current, the term 1% is the dissipative quasiparticle 
current which approaches the value V /R in the normal state, 
and 73 is the interference current. However, Eqs. d!6l > are more 
general, they include the Josephson oscillations of the spectral 
characteristics together with a nonequilibrium form of the dis- 
tribution function. 



IV. JUNCTION SPECTRUM 

To first order in eV, the spectral functions u(x,E,t) and 
v(x,E,t) obey static Usadel equations (II 31 with time-depen- 
dent boundary condition following from Eqs. ( Rai l and dl5al i. 



^0 (ud x v — vd x u) 



-4W(Mv) sin 2 O/2). 



(17) 



This boundary condition was found in Refs. [H] and[3 f° r a 
similar structure with the time-independent phase difference. 
Thus, at low enough voltages, the spectral characteristics of 
the junction adiabatically follow time variations of <p(t). 

Equations ( 1131 and dl7t can be supplemented by helpful 
identities following from the unity components of the matrix 
equations (j6ji and ( l8at . 

d t u = (iV/2)d x ((u,d x u) - (v,d x v)), (18a) 
$ ({u,d x u) - {v,d x v)) Q = -W(vg,cos0). (18b) 

In order to satisfy the condition ( U3bL we introduce usual 
parametrization u = cosh0, v = sinh0. Furthermore, we will 
neglect a small effect of suppression of the order parameter 
near the barrier, assuming A to be homogeneous [see comment 
to Eq. (12311. Then the equation and the boundary condition for 
the spectral angle 0, following from Eqs. (II 3ai and take 




FIG. 3: 2D profile of the DOS in the vicinity of the tunnel barrier at 
W = 0.05, at the moments of maximum depth of the Andreev band 
«t> = n). 



the form in a dimensionless variable z — x/^q, 



d-9 = 2A:sinh- 



k(E) = [iswh9 s (E)]- 



-27(f) sinh20 o , 7(0 = Wsin 2 [0(f) /2] 



(19) 
(20) 



The quantity 7(f) in Eq. ( 1201 has the meaning of a time- 
dependent depairing factor related to the discontinuity of the 
superconducting phase 0(f) at the tunnel barrier. When the 
phase approaches multiple of 2k, this factor turns to zero, and 
the spectral angle becomes homogeneous and equal to its bulk 
value 0s(E) = arctanh(A//i). At arbitrary (j), equations H9i 
and ( 1201 describe deviation of the spectral angle from 9s at 
the distance x ~ §0 from the barrier; thus, in a long junction, 
L <^o, we can apply the solution for a semi-infinite lead, 20 

tanh[(0 - 9 s )/4] = tanh[(0 o - B s )/4] exp(fe). (21) 

The equation for the boundary value 6o(E,t) of the spectral 
angle follows from Eqs. dl 9l - d2 11 . 



sinh 



9 s (E)-9 (E,t) _ 7(f) 
2 k{E) 



sinh 2 O (E,t). (22) 



The behavior of the DOS calculated by numerical solution 
of Eq. ( 1221 is shown in Figs. |2] and [3] At (j) = (we con- 
sider the phase within the period < <p < 2%), the DOS ap- 
proaches the BCS energy dependence depicted by the dashed 
line in Fig.|2ja). The current through the junction affects the 
singularity of DOS at the bulk energy gap, which becomes 
a beak-shaped peak with root singularities of the derivatives 
resulting from the divergence of the decay length k~ l (E) at 
E — > A [see Eq. (1191 1 . This divergence 20 is an analog of the 
long-range proximity effect at E — * in NS structures (see, 
e.g., a review 2 ^). The DOS identically turns to zero inside an 
oscillating minigap, at E < E*(t) < A. Similar conclusions 
can be drawn regarding behavior of the current spectral den- 
sities shown in Fig.|3a-c). 
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Within the subgap region, £*(f) < E < A, the DOS de- 
creases at the distance > E,q from the barrier (more precisely, 
at S^\E/A — 1|~ 1//4 , due to the long-range proximity effect 
mentioned above), as shown in Fig.|2jb). This implies that the 
subgap states form an oscillating cluster of bound states, "An- 
dreev band", which has triangular shape in the (£',x)-space 
as shown in Fig.[3](similar cluster appears near the lower gap 
edge, E = —A). The energy depth of the cluster is 
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^6- 



SE = A - w 6AW 4/5 sin 8 / 5 (0/2), 



(23) 



[see Eq. J28t > below], and it spreads over the distance of sev- 
eral |o from the tunnel barrier. We note that because of com- 
paratively large value of 8E, we may neglect in the spectral 
characteristics a smaller (~ AW) effect of local suppression 
of the order parameter A in the vicinity of the barrier. 12 Obvi- 
ously, the quantity 8E plays the role of the smallest character- 
istic energy of the spectrum in Eq. d!2i . 

The small value of the parameter W enables us to apply a 
perturbative approach for solving Eq. J22i . At the energies far 
enough from the gap edge, where 9 is of the order of unity, the 
quantity 9o is close to 9s, and Eq. 122\ leads to the following 
asymptotic relation for 9q, 



cosh O = cosh 9 S - 2 7(f) sinh20 s (/sinh 3 9s) 



1/2 



(24) 



However, when the energy approaches A, this expansion 
fails due to divergence of 9s at the bulk gap edge. This 
requires modification of the perturbative theory within the 
region \E — A\ <C Ar 2 - where the quantities 0o and 9s are 
large, while their difference, 9q — 9s, may have arbitrary 
value. Using these arguments, we hold only large expo- 
nents, exp 9o,s, in the hyperbolic functions in the right-hand 
side (rhs) of Eq. &22\ and use the asymptotic expression, 
exp 9s ~ \/2A/ (E — A) at small \E — A|. Then, introducing 
a dimensionless energy variable e and the normalized spectral 
function y(e) according to relations 

£ = p 2 (t)(E - A)/2A, p(t) = ^(f)] V5 , (25a) 
expOo =p(t)y(e), exp9 s = p{t)e- l l 2 , (25b) 

we reduce Eq. \22\ to a numerical algebraic equation fory(e), 

>2 



(26) 



According to Eqs. d25i . the parameter p(t) ~ 8E~ l l 2 deter- 
mines a characteristic scale of the spectral functions uq,vq w 
(1 /2) exp 0o m the vicinity of the bulk gap edge. 

The choice of the relevant complex root y(e) of Eq. i26\ is 
determined by the requirement for the asymptotic behavior at 
£ 3> 1 , which must coincide with the energy dependence given 
by direct perturbative expansion (I24> at 8E <C \E — A\ <C A. 
Within the main approximation, the function y(e) turns to 
zero at large £ as e -1 / 2 . At the energies £ smaller than 
£* = -(25/6) (2/3) */ 5 « -3.84, the function y(e) becomes 
imaginary, and therefore the spectral functions 
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FIG. 4: Spectral characteristic of the junction at W = 0.05, for dif- 
ferent values of the superconducting phase: = (curves indicated 
by 0), tt/4 (1), n/2 (2), and n (3). (a-c) - spectral densities l s , M 2 
and N 2 of the Josephson, interference and quasiparticle currents, re- 
spectively; (d) - spatially averaged DOS of the bound states. 

turn to zero at E < E*(t), where the minigap edge £*(f) is 
determined by relation following from Eq. (I25ai . 

£*(r)=A[l-C7 4/5 (?)], C = 25/(3 -6 1/5 )« 5.82. (28) 

We note that due to moderately small value W = 0.05 of 
the transparency parameter, the minigap values, obtained by 
numerical solution of Eq. i22\ and shown in Figs. [2] and 
13 slightly differ from their approximate values found from 
Eq. J28t > (by some 10%). Using the IPT approximation, how- 
ever, significantly simplifies both the analytical and numeri- 
cal calculations, allowing us to factorize the time and energy 
dependencies of the spectral functions. Applying this approx- 
imation to Eq. i2l\ . we obtain the spatial dependence of the 
spectral angle, 

exp9{x,E,t) = p(t)e- l/2 tanh 2 [s(e) - k(e,t)z/2], (29) 

tanh 2 s(£) =y(e)Ve, k(e,t) = [l^fe/ip{t)) 1/2 . (30) 

It is interesting to note that the spectral functions are inho- 
mogeneous within the leads at E = A and finite at any finite 
distance from the barrier, as it is seen in Fig. [3] which seems 
to contradict the divergence of the characteristic decay length 
k~ 1 (E) of the spectral angle at the gap edge. The explanation 
to this effect is the following: At E — > A (£ — > 0), the functions 
i(e) and k(e,t) in Eq. J30l > turn to zero as fi 1 / 4 , which cancels 
the divergence of the prefactor in Eq. d29l and results in the 
following expression for the spectral angle at the gap edge, 



N =M = ^p{t)Rey, I s = 



i/? 2 (f)Fmy 2 sin</>, (27) 



exp0(^A,f)=p(f)y(O)[l-z/v/2;>(f)y(O)] 2 . (31) 
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According to Eq. (13 11 . the spectral functions at E = A and 
z — > °o diverge, which restores their BCS-divergence at the 
gap edge within the bulk superconductor. 

V. DISTRIBUTION FUNCTION 

To first order in eV, equation @ in the Wigner representa- 
tion has the form of a diffusive kinetic equation^ 

Nd,f = Vd x (D+d x f), (32) 

where D + = (1/2) (1 + |m| 2 — |v| 2 ) is the dimensionless dif- 
fusion coefficient. At the reservoir, the distribution function 
approaches equilibrium value, f(—L) = n eq (E). The bound- 
ary condition at the tunnel barrier following from Eq. (I8b> . 

(D+d x f) =pd E fo, (33) 

determines the boundary value of the flow D+d x f of nonequi- 
librium quasiparticles escaping from the barrier. In this equa- 
tion, the function 

P{E,t)=eVI s (E,t)W/$ Q (34) 

describes the source of the nonequilibrium - the energy ex- 
change between the quasiparticles and the superfluid conden- 
sate due to the oscillations the supercurrent spectral density I s , 
This results in the DOS oscillations which can be expressed 
through I s by integrating the identity i ll 8 al l over the left lead, 
taking into account Eq. dl8bt and extending the integration 
over the whole negative semi-axis due to rapid decay of the 
time derivatives of all quantities at a small distance x ~ E,q, 

,o 

(d,N) x (E,t) = J d,N(x,E,t)dx = 'DdEP(E,t). (35) 

Due to the existence of oscillating Andreev band in the 
vicinity of the barrier, the behavior of the quasiparticles 
strongly depends on their energy. Indeed, most of excitations 
with the energy E > A spend a short time ~ A~ 1 near the bar- 
rier, they rapidly diffuse away along the leads and escape into 
the reservoirs (travelling quasiparticles). The quasiparticles 
with low energy approaching the contact at the distance ~ <^o 
while the depth of the Andreev band increases (0 < <j> < 7C) 
are trapped by the Andreev band and spend much larger time 
~ (eV) -1 inside the band, following the spectrum oscilla- 
tions. During the next half-period (n < <j> < 2%), the depth 
of the Andreev band decreases, and the trapped quasiparticles 
are pushed out to the extended states. Such a physical pic- 
ture is similar to the quasiparticle dynamics within the surface 
skin layer of a superconductor irradiated by rf electromagnetic 
field. 22 Below we will perform separate analysis of the kinet- 
ics of the travelling and trapped quasiparticles. 

A. Travelling quasiparticles 

At E > A, the coefficients N and D + in Eq. 1321 rapidly 
vary in the near vicinity x ~ of the barrier and then, 



at i > approach their values in a bulk superconductor, 
N -» N S (E) = E/\/E 2 -A 2 , D+ -> 1. Correspondingly, the 
solution of Eq. 1321 has the form of a slowly varying func- 
tion /(°) with a small but rapidly varying addition /W, which 
vanishes at the distances i>^o from the barrier. Thus, within 
the main approximation, the population of the extended states, 
/> = f(°\x,E > A,f), satisfies Eq. (I32t with the asymptotic 
values of the coefficients at x 3> 

^«9,/ > = D«9 2 / > . (36) 

To derive effective boundary condition to this equation, we 
subtract Eq. (1361 from the initial kinetic equation (1321 . keep- 
ing large spatial derivative of the function 

(N - N S ) d t h = Dd x [(D+ - 1 )d x f > + D+d x fW] , (37) 

and then integrate Eq. (I37l > over x along the left lead. Since 
all terms in this equation vanish at x ^ §q, we extend the in- 
tegration over the whole semi-axis, similar to Eq. d35l >. and 
take the smooth function /> in its left-hand side (lhs) at x = 0. 
Then, using Eq. (I33> . we obtain the boundary condition, 

(N-N s ) x (d,f>)o = D (j3«9 £ /> - d,/>) . (38) 

The distribution function /> consists of a static part g(x,E) 
and a dynamic (time-dependent) part, />(x,£,r) = g(x,E) + 
h(x,E,t). The static part linearly varies along the lead and ap- 
proaches equilibrium distribution function n eq (E) at x = —L, 

g(x,E)=go(E) + (x/L)[g (E)-n eg (E)}. (39) 
The dynamic component h(x,E,t) has the form, 

h( x ,E,t) - Y,,M E ) e ~ ineVt+K " x i W' = °> < 4 °) 

K n (E) = y / |n|ey/Y 5 (£)2)- 1 exp[-(/7i:/4)sgnn]. 

Here (. . .) t denotes time averaging over the period T = 2n/eV, 
and h n (E) are the Fourier harmonics of the function h(Q,E,t), 

r% At 

h n {E) = / —h(Q,E,t)e meVt . 
Jo X 

Equation J40i was obtained assuming the decay length Ly — 
K„ ~ y/D/eV of the oscillations of the distribution function 
to be much smaller than the junction length L; this enables 
us to apply the solution for a semi-infinite lead. In the op- 
posite case, Ly S> L, the reservoirs effectively suppress these 
oscillations, and the function h rapidly decays while the volt- 
age decreases. The equation for the harmonics h n (E) follows 
from Eq. d38l. 

T>- l (N-N s ) x d t h{0,E 7 t)-pd E h(0,E,t) + {d xg ) 

+ Y jn h„(E)K„e- ineVt = J3 d Eg0 (E). (41) 

The first two terms in the lhs of Eq. J41I . proportional to small 
applied voltage eV , can be neglected as compared to the fourth 



7 



term. Then, averaging Eq. (1411 over t and using the identity 
( I35i . we obtain the expression for the third term, 

{d x g) = d E {Ph)u (42) 

which therefore is also proportional to small eV and can be 
neglected as well. As the result, the approximate solution 
of Eq. J41t reads h n = fi„K~ 1 d E go, i.e., the dynamic part 
h(x,E,t) of the distribution function is expressed through the 
boundary value of the static part, go(E). 

It follows from this consideration that the dynamic part of 
the distribution is generated by the oscillations of the quasi- 
particle spectrum at E > A. These oscillations provide the 
energy transfer from the electric field to quasiparticles, which 
generally results in the 'heating' of quasiparticles, i.e., their 
redistribution towards higher energies (pumping effect). The 
heating is described by a diffusion equation for go in the en- 
ergy space following from the Eqs. (142 \ and (I39L 

d E {D E d E go) = L(d x g) = go~ n eq , (43) 
D E =L-£ n \p n \ 2 K-\ 

Similar equation has been obtained in Refs. |2J for MAR-in- 
duced heating in long SNS junctions, where the nonequilib- 
rium is constrained by inelastic collisions. In our strongly in- 
homogeneous case, the role of relaxation factor in Eq. J43i is 
played by the diffusive flow (d x g)o of nonequilibrium quasi- 
particles from the junction. The intensity of the heating ef- 
fect is determined by the estimate of the diffusion term in 
Eq. (I43> with the equilibrium function n eq , d E (D E d E n eq ) ~ 
(L/£o)W(eV /8E) 3 / 2 . Due to the presence of the two small 
factors, this term is small for reasonable lengths of the leads, 
which enables us to approximate the function go with the equi- 
librium distribution function, go « n eq . Then the boundary 
value of the distribution function of travelling quasiparticles 
reads, 

f > (0,E,t)=n eq (E)+h(0,E,t), (44) 
h(0,E,t)=n' eq 1 £ n PnK- 1 e- ineVt , n' eq = d E n eq {E). (45) 

B. Trapped quasiparticles 

At small voltages, eV <C A, the spatial size £,o of the An- 
dreev band is much smaller than the smallest kinetic length 
Ly, therefore the main part /(°) of the distribution function of 
the trapped quasiparticles is spatially uniform, /< « /w [E < 
A,t). Then the kinetic equation d32t takes the form, 

Nd t f < = Vd x (D+djM). (46) 

By averaging Eq. J46i over x and using the boundary condi- 
tion J33I . we obtain partial differential equation (PDE) for the 
function /< , 

(N)Af<-Vpd E f < =0. (47) 

As long as the coefficients of this PDE satisfy the identity d35l . 
the equation d£, = (N) x dE + CD/3 dt = for its characteristic 
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FIG. 5: Levels of equal number of states £,{E,t) = const vs phase 
difference at W = 0.05. Bold line depicts the edge £»(/) = E{t, = 
0,t) of the minigap, and to indicates a particular moment of trapping 
to the trajectory § = 0.8A. 

curves £(E,t) — const in the (E,t) space is easily integrated, 

= V 1 / dE'(N(x,E',t)) x . (48) 

JE t (f) 

This allows us to reduce the PDE d47l to the ordinary differ- 
ential equation for the function /< along the characteristics, 

d,f<(^t)\^ comt = 0. (49) 

The quantity %(E,t) has the meaning of averaged number 
of states with the energy smaller than given energy E, nor- 
malized ov er N E ; in the bulk superconductor, it approaches 
®(E - A)V£ 2 -A 2 . The bottom E*(t) of the Andreev band 
represents the reference point, £, [E* (t),t] =0. Within the IPT 
approximation, Eq. J48b reads 

|(£,f) = A[47(f)] 1/5 7(e), 7(e) = J* de'n(e'), (50a) 

n(e) = |er 1 / 2 l m y^(i), (50b) 

where the function n(e) is related to the averaged DOS as 
(N) z = [2y i {t)]- l l 5 n{e). As follows from Eq. (BObt . the func- 
tion (N) Z (E) plotted in Fig.|4ld) exhibits a singularity at the 
bulk gap edge, due to the long-range proximity effect in the 
superconductor [see comments to Eq. (I22H . 

According to Eq. ( I49K the distribution of trapped excita- 
tions is a function of the "integral of motion" t, (E,t), 

f<(E,t)=f[S(E,t)], (51) 

and therefore holds a constant value along the trajectories 
£,(E,t) = const in the (E,t) space shown in Fig. EI These tra- 
jectories can be interpreted as diffusive analogue of Andreev 
levels adiabatically moving in a slowly varying external field 
(cf. Ref. 1241) . From this we conclude that the trapped quasi- 
particles are completely dragged by the spectrum oscillations, 
and their energy periodically changes, in accordance with the 
relation following from Eqs. (I48i and ( 1351 . 

d t E\ i=aBA = -'Dpiff)-\ (52) 
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Since the bound states at the tunnel barrier appear and van- 
ish periodically, their population is imposed by the distribu- 
tion of travelling quasiparticles at the bulk gap edge. This is 
expressed by the matching condition 



/<(§)=/>[0,A,<b(S)], 0<(j><7C. 



(53) 



Here ?o(4) is the moment of the crossing the gap edge by 
quasiparticle with the given value of | (see Fig. |5Jl, i.e., it 
is the smallest solution of the equation E(% ,fo) = A within the 
period of the DOS oscillations. 

Since the trapped quasiparticles spend long time at the bar- 
rier, the heating effect for them becomes well pronounced, in 
contrast to the travelling quasiparticles. Indeed, as follows 
from Eq. (I45L the time-dependent part h of the distribution 
function turns to zero at E = A due to singularity of Ns, and 
therefore the distribution function of trapped quasiparticles 
turns into a plateau (i.e., does not depend on energy), 



/<(§)=So(A) 



,(A). 



(54) 



Non-adiabatic correction to the distribution function d54t 
is produced by inelastic relaxation. To include the effect of 
inelastic collisions on the quasiparticle kinetics, we will add 
a collision term in the T-approximation, Nt^ 1 [n eq — /], to the 
rhs of the kinetic equation \A6\ . Within such a model, the 
kinetic equation for the trapped quasiparticles reads 

d,U = t- 1 {n eq [E{t; , t )] - /< (4 , t ) } . (55) 

The initial condition to Eq. fl55i follows from Eq. j53> . 



MM 



(A). 



(56) 



VI. ELECTRIC CURRENT 

Proceeding to the calculation of the electric current dl6at . 
we start from the analysis of the ac Josephson current I\ sin0. 
First, we note that the effect of the DOS oscillations on 
the amplitude I\ , calculated with the equilibrium distribution 
function n eq , gives rise only to small corrections to the tunnel 
model result, l[ q = {KA/2eR)\.<m\\(A/2T)& Indeed, in this 
case, the integral in Eq. (I16bt is determined by the imaginary 
Matsubara energies ico„ = KiT (2n + 1 ) far from the gap edges, 
which allows us to apply ordinary perturbation approach [see 
Eq. J24H to the calculation of the spectral density I s . This re- 
sults in the correction of the order of W to Zf ? j 12 which will be 
neglected below. 

In what follows, we will focus on the nontrivial contribu- 
tions of nonequilibrium quasiparticles to the Josephson cur- 
rent Zj. Such contributions, Z/,(f) and Z<(f), come from the 
dynamic part h of the distribution of travelling quasiparticles 
and the distribution /< of trapped quasiparticles, 

'*(') = \ E e ~'" eVt I" dEI s K- l ^ n n' eq {E), (57a) 
r A 



Due to rapid convergence of the integral in Eq. ( I57al > at £ — 
A < 8E <C (A, T), the function n' (E) can be taken at E = A. 
For similar reason, we apply the IPT expression d27t for the 
spectral density I s , which leads to the following result, 



where 

Fi(T) 
P(t) 

Ci 



w eR V 2A w r ' 



, , . A A 

:A «e 9 ( A ) = ^ COsh 

d sin 2v - £ {-l) m Vmsm{m<l> - n/4) 



(58) 



2 ^ (m 2 - v 2 ) B(v + m, v - m) ' 

: 2 : - v / 2 rrfee 1 / 4 [Im 3 ; 2 (e)] 2 « 0.15. v = 0.2, 
Jo 



and B(x,y) is the Euler's beta function. The time average of 
the current (15 8> does not vanish, 



C 3 = (P(t) sin (j>), =0.018, 



(59) 



and, moreover, it exhibits a strongly nonlinear (~ ^/V) voltage 
dependence. We recall that the validity of these results is re- 
stricted by the condition of small diffusion length as compared 
to the junction length, L v = ^D/eV <C L, or, equivalently, 
eV Efh, where Ejj, = D/L 2 is the Thouless energy [see 
comments to Eq. J40H . In the opposite case, eV £Y/,, the 
oscillations of the distribution function are damped by equilib- 
rium reservoirs, and the nonequilibrium dc current, produced 
by travelling quasiparticles, rapidly vanishes. Similar damp- 
ing effect is caused by inelastic collisions, when the inelastic 
scattering length becomes smaller than the diffusion length 
L V . 

The dissipative dc current of travelling quasiparticles re- 
sults from non-adiabatic (diffusive) evolution of the distri- 
bution function h(t). Similar conclusion is also true for the 
trapped quasiparticles: only the non-adiabatic part of the dis- 
tribution function contributes to the dissipative dc current. 
Indeed, inserting the adiabatic distribution function n eq (A), 
Eq. J54> . into Eq. J57bi . and taking advantage of a small inte- 
gration interval, which allows us to expand the difference of 
the equilibrium functions over E — A, we get, 



&( t ) = ^ic 7 W 4 / 5 sin 8 / 5 ^sin</>, (60) 
< w eR - 2 r 







C 2 



>-2V 



dee\my 2 {e) «3.11. 



This equation contains only odd harmonics of the ac Joseph- 
son current. The dc current results from the non-adiabatic cor- 
rection 8f to the distribution function, found from the solution 
of the kinetic equation j55l >. 



If 1 / f \ AFi 

I<(t) = ^ R -J E dE W <[!(£, t)A - n eq (E)}. (57b) J* = — (y dEI s 8fg(E,t),t})= -^-W^ 5 K(eV). (61) 
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The factor K(eV) in Eq. (16 1 i is a complicated numerical func- 
tion of the relaxation parameter eVz £ . In the limits of weak 
(eVT e ^> 1) and strong {eVi £ C 1) relaxation, this function 
decreases as (eVT £ )~ l and eVx £ , respectively. In the spirit 
of our modelling approach to the problem of inelastic scatter- 
ing, it is reasonable to approximate K(eV) by usual relaxation 
factor 



K(eV) : 



aieVx £ 
a 2 + {eVr £ 



(62) 



where the parameters a,\ 2 are to be evaluated from the asymp- 
totics of K(eV) in both limiting cases. 

In the weak relaxation limit, eVi £ ^> 1, the adiabatic dis- 
tribution function is given by Eq. (|54}, ff = n eq (A); then the 
non-adiabatic correction 8f is determined by the local energy 
E(^,t) averaged along the trajectories £, = const, 

8f£,t) = z- 1 n' (A) f dt'[E(^t')-A\. 



Substituting 8f into the expression (00) and using the IPT re- 
lations ( 15 Oi l for the function % (E,t), we obtain the coefficient 
a,\ in the relaxation factor K(eV), 

a, = B (°-y- 3 ) [°delmy 2 (e) f 1 dxeo\x^ 5 J{e)] « 1.38, 

2 z / ;> ^ A* Jo 

where £o(C) is the solution of equation £ = J(so), and the 
function 7(e) is defined in Eq. ( I50a> . 

In the opposite limit of strong relaxation, eVx £ -C 1, the 
initial condition d56i is quickly 'forgotten', and the adiabatic 
part of the distribution is the equilibrium function of the lo- 
cal energy, = n eq [E(E,,t)\. In this case, the non-adiabatic 
addition 8f is proportional to the time derivative of E(£,,t) 
along the trajectory, 

8f(^t)^-T £ n' eq (A)d t E(^t). 

Using the expression i52\ for d t E and calculating I< in 
Eq. (16 It in the IPT approximation, we arrive at the following 
relation between the coefficients a,\ and a 2 , 



<Xi _ f 71 afTsin z 2T 
« 2 Jo 2 7 /%sin 2 / 5 



de 



[Imy 2 (e)] 2 



1.74. 



Now we turn to evaluation of the quasiparticle current I 2 
and the interference current I^cos^). These currents are pro- 
portional to small applied voltage, and therefore they can be 
calculated within the equilibrium approximation for the quasi- 
particle distribution neglecting small nonequilibrium correc- 
tion. The contribution of the trapped quasiparticles is small, 
because of the small occupied phase volume, and can be ne- 
glected. Moreover, in the limit of weak relaxation, this con- 
tribution identically turns to zero because of energy-indepen- 
dent distribution /< . Within the BCS approximation for the 
spectral functions Nq and Mo, the integrals in Eqs. J16cl > and 
d!6dt logarithmically diverge at the gap edge. The effect of the 
phase difference on the quasiparticle spectrum eliminates this 
divergence and effectively cuts the spectral functions at the 



value {8E)~ X I 2 [see comments to Eq. ( I25H . which is equiva- 
lent to the effective cut-off in energy, E — A > 8E. Using the 
IPT expressions d27l > for the spectral functions at E — A <C A, 
we obtain to the main order in W, 



h(t)=h(t) + j[l-n eq (A)], 



h(t) = j[Fi\n S fa/f/5(t)-F 2 



(63) 
(64) 



where 



A 2 dE 



E 2_ A 2 K,(^)-n' eq (E)]; 



In a— de 

Jo 



2 ^ ©(£-!) 



Re 2 y(e) 



2 

-ln2; 



0.13. 



We note that within the tunnel model* 2 ^ the role of the cut- 
off factor is played by the quantity eV, which enters large 
logarithmic factor in the expression for the current. In our 
case, according to the adiabaticity criterion H2\ . the effec- 
tive cut-off factor, y 4 / 5 ~ SE, is much larger. From this we 
conclude that within the region of applicability of the adia- 
batic approach, eV <C 8E, the quasiparticle and the interfer- 
ence currents are logarithmically smaller as compared to the 
results of the tunnel model. While the voltage increases and 
exceeds the depth of the Andreev band, eV > SE, the tunnel 
model approximation for the currents I 2 and It, becomes valid. 

Due to the presence of a time-dependent factor y(f ) in the 
logarithmic terms, the currents I 2 and It, oscillate in time and 
exhibit logarithmic singularities when 7(f) turns to zero, i.e., 
at <p = 0. In the vicinity of these points, the Andreev band van- 
ishes, which violates the condition of adiabaticity i ll 2t . Fol- 
lowing remarks to the equations d63l > and i64l . we can semi- 
quantitatively describe the whole shape of the oscillations of 
the quasiparticle and interference currents by adding a small 
regularization term eV /A to the denominator of the logarith- 
mic argument in Eq. J64l >. 

By averaging Eq. d63l over time, we obtain the dissipative 
part of the current I 2 (t), 



(h) t = ^[F l ln^b/W 4 / 5 +F 3 



\nb = \na 



7Z 



lnshT 8/5 T= 1.24, 



F 3 (T) = l-n eq (A)-F 2 (T). 

The interference current also has a dc component due to the 
oscillations in ^3, though its magnitude, 

V 

(I 3 cos <j)) t = — FilnVc, 
K 

i*7l At 

]nc= / — cos2Tlnshr 8/5 T = 0.8, 
Jo n 

is small with respect to the logarithmic term in (I 2 ). The sum 
of the quasiparticle and interference dc currents is given by 
expression, 

4l = (I 2 + h cos </>), = £ [Fi ln(2.77/W 2 / 5 ) + F 3 ] . (65) 

K 
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FIG. 6: Temperature dependencies of the functions F\, Fi, and F3. 

The temperature dependencies of the functions F\(T), F^{T) 
and Fj(T) are shown in Fig. [6] 

VII. DISCUSSION 

Analyzing different calculated current components, we 
conclude that the most significant deviation from the results 
of the tunnel model comes from the nonequilibrium Josephson 
current, while the quasiparticle and interference currents un- 
dergo minor (logarithmic) changes. There are the two features 
to be mentioned. The first is a non-monotonous voltage de- 
pendence of the current of trapped quasiparticles described by 
the relaxation factor K(eV) [see Eqs. i61\ and J62H . which has 
a maximum at eV ~ TjT 1 . Due to rather large value of the pa- 
rameter t £ A > 10 2 in most superconductors, this current may 
exceed linear contributions of the quasiparticle and interfer- 
ence currents at moderately small values of the transparency 
parameter W. In this case, the I-V characteristic exhibits an N- 
like feature, as shown in Fig. 0a), and the linear conductance 
at small bias considerably (by the factor W 4 ^ 5 T £ A) exceeds the 
tunnel model conductance, which is relevant for larger bias 
[zero bias conductance peak, Fig. 0b)]. 

The second feature is a nonlinear, proportional to y/V, con- 
tribution to the current-voltage characteristics produced by 
the travelling quasiparticles. We note, however, that because 
of small numerical value of the constant C3 in Eq. d59l . the 
crossover of I(V) to nonlinear behavior for reasonable val- 
ues of W actually occurs at very small voltage. The reason 
for this is a numerically small magnitude of the supercurrent 
spectral density I S (E) above the bulk gap edge, as is obvi- 
ous from Fig.|4la). Physically, this is due to a rapid decrease 
of the probability of the "over-the-barrier" Andreev reflection 
(reflection at the energy outside the energy gap). As the result, 
the density of nonequilibrium quasiparticles produced by the 
oscillations of I s at E > A appears to be small as compared 
with the "natural" width 8E of the perturbed spectral region, 
which results in smaller contribution of the travelling quasi- 
particles compared to the trapped ones. 

In order to estimate the characteristic parameters of the 
junction for the transparency parameter W = 0.05 accepted in 
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a) b) 

FIG. 7: dc current (a) and differential conductance (b) vs voltage for 
different values of T e A = 50 (1), 100 (2), and 200 (3). All curves are 
plotted at W = 0.05, T = 0.5T C . 

our numerical calculations, we will assume the junction area 
to be 200 x 200 nm and the thickness of the leads as well as 
the elastic scattering length to be 50 nm. For Al leads, this 
results in the sheet resistance /?□ ss 0.3 fi and ^(^0) ~ 0.45 
Q. at ^0 ~ 300 nm. Then, according to Eq. l|9}, the tunnel- 
ing probability, the junction resistance and the critical current 
approach the values T w 0.01, R « 4.5 £2, and I c w 70 /J.A, 
respectively. Thus, the characteristic voltage region, at which 
the nonequilibrium dc current dominates, is of the order of 
several microvolts. 

It follows from the presented estimates and also from Fig.0 
that the nonequilibrium effect appears at rather small applied 
voltage, which makes it difficult to observe in practice because 
of the jumps to the Josephson branch. To facilitate the ob- 
servation, the net Josephson current must be suppressed by 
applying magnetic field or using the dc SQUID setup. Re- 
markably, the effect survives even in the absence of the net 
Josephson current: Local Josephson currents flowing through 
different junctions (or parts of the junction) generate dissi- 
pative dc current flowing in the same direction, which is de- 
termined by the applied voltage (i.e., time derivative of the 
phase) rather than the local values of the phase difference. In- 
deed, let us consider two junctions connected in parallel (dc 
SQUID) and apply half-integer magnetic flux. Then the phase 
differences at the junctions will be shifted by %. This will lead 
to different signs of the current spectral density in Eq. J61t for 
different junctions, however, the nonequilibrium distribution 
function will also change the sign, which can be proven by di- 
rect calculation. Actually, the equality of dc currents in both 
junctions follows from a simple fact that the constant phase 
shift is equivalent to the change of the time reference point, 
which obviously cannot affect the value of the time-averaged 
(dc) current. Thus, the net dc current will double while the 
ac current disappears. Such situation resembles the case of 
a long SNS junction, 10,24 where the Josephson effect is sup- 
pressed due to decay of the superconducting correlations in- 
side the normal part of the junction, and the nonequilibrium 
effect can be observed at very small applied voltage. 

It is instructive to compare the mechanism of nonequilib- 
rium tunnel current in diffusive junctions with mesoscopic 
picture of superconductive tunneling in point contacts given 
by the MAR theory. 5,7 ' 8 It has been already mentioned in the 
Introduction, that the Andreev band is the qualitative analog 
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of the Andreev bound levels in point contacts, the both are 
formed by the same physical mechanism - Andreev reflec- 
tion by the jump of the superconducting phase at the junc- 
tion. The time oscillations of the Andreev levels in volt- 
age biased point contacts, and the quasiparticle exchange be- 
tween the levels and the continuum is the adiabatic descrip- 
tion in the time domain of the coherent MAR. 8 The latter 
also include the quasiparticle transitions across the energy 
gap, which in the real-time picture correspond to the Landau- 
Zener tunneling between the Andreev levels (bands). 7 Such 
tunneling provides the flow of quasiparticles along the energy 
axis (spectral flow), proportional to the gradient of the dis- 
tribution function in the energy space, dgf, within the sub- 
gap region, \E\ < A. These processes are weak at small volt- 
age compared to the continuum-bound level exchange, and 
they are neglected within the adiabatic approximation adopted 
in the present paper. For this reason, the spectral flow of 
the trapped quasiparticles is blocked, which results in a fiat 
(energy-independent) distribution ( 1541 of these quasiparticles 
in the absence of inelastic relaxation. 

However, there is considerable quantitative difference be- 
tween the point contacts and diffusive junctions: In the point 
contacts, most of the quasiparticles reflected by the tunnel 
barrier escape to the reservoir, while in the long diffusive 
lead the quasiparticles multiply collide with the barrier due 
to the impurity backscattering. This essentially increases the 
probability of the coherent tunneling and Andreev reflection 25 
and, correspondingly, enhances the effect of the phase differ- 
ence on the junction spectrum. As the result, the depth of 
the Andreev band in our long-arm geometry, 8E ~ AW 4//5 in 
Eq. d23i . considerably exceeds the depth ~ Ar of the Andreev 
level in a point contact with comparable transparency. 

It is interesting to mention that generally the role of the 
junction transparency in tunnel junctions with diffusive elec- 
trodes is played by the parameter W rather than the barrier 
transparency Indeed, comparison of the right- and left- 

hand sides of Eqs. (|8} shows that as long as the 'natural' scale 
of the currents J and J is proportional to the gradients of the 
Green's and distribution functions in the vicinity of the bar- 
rier, ~ it is the magnitude of W (not T itself) which de- 
termines the effective barrier strength. For example, at large 
IV>1, the commutators in Eqs. (|8| are to be small, which 
implies continuity of all functions at the barrier. From this 
we conclude that at W 3> 1 , the phase and voltage are contin- 
uously distributed along the whole structure, and the barrier 
does not affect the current transport even if its transparency is 
small, r< 1, provided T 3> I /4o- In other words, at large W, 
the critical current of the tunnel junction, formally estimated 
by the Ambegaokar-Baratoff formula for a low-transparent 
junction, 26 exceeds the critical current of diffusive supercon- 
ducting leads. Physically, this effective 'blooming' of an 



opaque barrier results from the multiple coherent backscat- 
tering of quasiparticles by impurities mentioned above. 25 
VIII. CONCLUSIONS 

We have studied the current transport through the voltage 
biased tunnel junction with long diffusive superconducting 
leads at low applied voltage, eV <C 2A, and finite tempera- 
tures. In contrast to the tunnel model* 2 ^ we consider the non- 
equilibrium effects in the spectrum and distribution of quasi- 
particles in the vicinity of the barrier. Using the small value 
of the Josephson frequency with respect to other characteris- 
tic energies of the problem, we apply a quasi-adiabatic ap- 
proach to the analysis of the time evolution of quasiparti- 
cles. Within such approach, we obtain a physically trans- 
parent picture of the quasiparticle spectrum adiabatically fol- 
lowing the time-dependent difference of the superconducting 
phases. This results in local oscillations of the spectrum of 
travelling quasiparticles (E > A) and formation of oscillating 
Andreev bound states (Andreev band) within the subgap re- 
gion, E*(t) < E < A, at the distance of the order of the coher- 
ence length near the barrier. The quasiparticles trapped by 
the Andreev band are completely dragged by oscillations of 
the junction spectrum which reflects complete multiple An- 
dreev reflection (MAR) of the subgap quasiparticles and re- 
sults in generation of higher odd harmonics of the ac Joseph- 
son current. The inelastic relaxation of the trapped quasiparti- 
cles produces a non-adiabatic component of their distribution 
which manifests itself in the nonequilibrium contribution to 
the dc current. At low enough voltages, this contribution may 
considerably exceed the quasiparticle dc current given by the 
tunnel model; by this reason, the resulting I-V characteristic 
shows 3\f-like feature, with the maximum at eV ~ TjT 1 . The 
travelling quasiparticles also deviate from equilibrium due to 
partial drag (which corresponds to the over-the-barrier MAR) 
confined by their fast diffusion from the barrier. This results 
in the additional contribution to the dc current, proportional 
to y/V. We note that our approach can be easily extended to 
the current bias regime; in this case, arbitrary time-dependent 
phase must be assumed in the calculation. 

Effect of travelling quasiparticles is interesting, in principle 
it dominates at small voltage because of the non-analytical 
voltage dependence. Unfortunately, this contribution is nu- 
merically small in the case of the tunnel barrier considered 
here, because of the small spectral density of the Josephson 
current above the gap. However, there are no fundamental rea- 
sons to expect this contribution to be small in other junctions. 
It would be interesting to find favorable junction geometry. 
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